NASA TECHNICAL NOTE 




■ 






XfWL TECHNICAL LIfc, 

I< j r?Ti * V ^ . N . M . 


COMPARISON OF TECHNIQUES 
FOR APPROXIMATING OCEAN BOTTOM 
TOPOGRAPHY IN A WAVE-REFRACTION 
COMPUTER MODEL 


Lament R. Poole 

Langley Research Center 
Hampton, Va. 23665 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 



I ■ , ’ 

WASHIN6TON/d. C. • NOVEMBER 1975 



tech library kafb, nm 



0133^01 


1. Report No. 2. Government Accession No. 

NASA TN D-8050 

3. Recipient's Catalog No. 

4. Title and Subtitle 

COMPARISON OF TECHNIQUES FOR APPROXIMATING OCEAN 
BOTTOM TOPOGRAPHY IN A WAVE-REFRACTION COMPUTER 
MODEL 

5. Report Date 

November 1975 

6. Performing Organization Code 

7. Author(s) 

Lament R. Poole 

8. Performing Organization Report No. 

L-10320 

10. Work Unit No. 

161-03-03-01 

9. Performing Organization Name and Address 

NASA Langley Research Center 
Hampton, Va. 23665 

11. Contract or Grant No. 

13. Type of Report and Period Covered 

Technical Note 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 

14. Sponsoring Agency Code 

15. Supplementary Notes 

Appendix entitled “Description of Contouring Technique” by James F. Kibler. 

16. Abstract 


A study of the effects of using different methods for approximating bottom topography in 
a wave-refraction computer model was conducted. Approximation techniques involving quadratic 
least squares, cubic least squares, and constrained bicubic polynomial interpolation were compared 
for computed wave patterns and parameters in the region of Saco Bay, Maine. Although substantial 
local differences can be attributed to use of the different approximation techniques, results indi- 
cated that overall computed wave patterns and parameter distributions were quite similar. More 
extensive quantitative comparison of computed results with actual wave data would be required in 
order to select the best approximation technique for refraction studies in which local accuracy in 
computed results is critical. 


17. Key Words (Suggested by Author(s) ) 

Water waves 
Wave refraction 
Numerical approximation 


18. Distribution Statement 

Unclassified — Unlimited 


Subject Category 48 


19. Security Classif. (of this report) 

Unclassified 


20. Security Classif. (of this page) 

Unclassified 


21. No. of Pages 

44 


22. Price* 

$3.75 


For sale by the National Technical Information Service, Springfield, Virginia 22161 




COMPARISON OF TECHNIQUES FOR APPROXIMATING OCEAN BOTTOM 
TOPOGRAPHY IN A WAVE-REFRACTION COMPUTER MODEL 


Lament R. Poole 
Langley Research Center 

SUMMARY 

A study of the effects of using different techniques for approximating bottom 
topography in a wave-refraction computer model was conducted. A mid-Atlantic 
continental-shelf refraction program was modified to study the region of Saco Bay, Maine. 
The program employed, on option, topography approximation techniques involving a 
quadratic least-squares approach, a cubic least-squares approach, or constrained bicubic 
polynomial interpolation. The results were computed using each of the techniques in 
conjunction with a given set of initial wave conditions. Comparisons of the results indi- 
cated that overall wave-ray and crest patterns were not strongly affected by the choice of 
approximation technique. The computed crest patterns compared equally well with the 
pattern observed in vertical aerial photographs of a Saco Bay wavefield. Differences among 
the comparative patterns in the paths of corresponding individual rays were noted. The 
path of one particular ray, which was computed using the bicubic polynomial interpolation 
technique, changed abruptly as the ray passed near an isolated extremum in the input 
bathymetry data field. Comparative contours at selected levels of wavelength were quite 
similar; comparative contours at selected wave-height levels, while exhibiting considerable 
local differences, also showed similarity in overall orientation. For applications in which 
local accuracy in computed results is critical, more extensive comparison of computed 
results with actual data on wavelength, direction, and height would be required in order to 
select the best technique for topography approximation. 

INTRODUCTION 

In the past, man has made extensive use of the world’s coastal regions and adjacent 
continental-shelf waters for commerce, recreation, and defense purposes. In view of 
proposals for exploitation of offshore petroleum deposits, construction of offshore power 
and port facilities, and offshore dumping of sewage sludge and other wastes, the future 
can only lead to more extensive and complex interactions with the coastal marine 
environment. Identification of potential safety and pollution problems arising from these 


interactions requires improved experimental and analytical techniques for studying the 
physical phenomena involved. 

One such physical phenomenon, a knowledge of which is a most important element 
in the processes of both offshore and shoreline facility design and site selection, is the 
behavior of surface waves as they propagate over the continental shelf and impinge on the 
shoreline. The analytical tool necessary for the study of this process is the wave-refraction 
model, one such as those models described in references 1 and 2. In the past, such models 
have been generally limited to small geographical areas such as harbor entrances. Recently, 
however, the Langley Research Center and the Virginia Institute of Marine Science, in a 
cooperative effort, have developed a computerized model (ref. 3) for study of wave 
refraction on the mid-Atlantic continental shelf from Cape Hatteras to Cape Henlopen. 

A model of this scale (100 n. mi. by 210 n. mi.) is needed in order to identify regional 
differences in wave behavior to aid in site selection for offshore and shoreline facilities. 

The variables of primary influence in an ocean wave-refraction study are the local 
water depth and its spatial variation. The bottom topography for any large-scale 
geographic region varies in a generally random fashion. Therefore, the only practical 
method for the input of bathymetry data to a refraction model is through some type of 
spatial grid format in which known depth values are supplied at the nodes of the grid. 

In turn, use of such a grid-type input for bathymetry data necessitates some procedure for 
determining depth values at intermediate points between nodes and for computing spatial 
derivatives of the depth for input to other refraction calculations. 

Two procedures which have been widely used to approximate bottom topography in 
refraction computer programs are the quadratic least-squares technique (refs. 1 and 3) and 
the cubic least-squares technique (ref. 2). Another approximation technique, which has 
been used by the Coastal Engineering Laboratory of the United States Army Corps of 
Engineers, involves constrained interpolation with a bicubic (ref. 4, p. 97) polynomial.^ 

Each of these techniques is easily applied in a mosaic fashion to large-scale geographic 
regions. While the constrained bicubic interpolation ensures global continuity in the depth 
values and certain partial derivatives, the least-squares procedures do not insure continuity 
in either the depth or its derivatives. Such continuity is desirable from the standpoint of 
uniqueness in the solution of the refraction equations. On the other hand, the least-squares 
techniques do provide some smoothing of input bathymetry data (thereby possibly averaging 
out random error in the data), whereas the constrained bicubic polynomial interpolates 
through the input data without such smoothing. 


^This interpolation algorithm was obtained through private communication from 
Dr. D. Lee Harris, Coastal Engineering Research Center, Fort Belvoir, Virginia. 
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Because of these inherent differences, a study was conducted to investigate the effects 
on wave-refraction calculations arising from the use of the different techniques to approximate 
bottom topography. A modified version of the computerized wave-refraction model of 
reference 3, which employs the quadratic least-squares approximation method, was used in the 
study. The cubic least-squares and constrained bicubic polynomial approximation techniques 
were also included in the model. In addition, bathymetry data for the region of Saco Bay, 
Maine (provided by Dr. Stewart C. Farrell, Stockton State College, Pomona, New Jersey) were 
substituted for the mid-Atlantic continental- shelf data in the basic model. This small 
geographic region was selected in order to reference computed wave patterns to vertical aerial 
photographs of an actual Saco Bay wave pattern (ref. 3, p. 11), as well as to facilitate 
comparison of overall patterns computed using the different approximation techniques. 

A brief discussion is given of the bottom-topography approximation techniques as 
applied mosaically to a fixed grid system. For a set of initial wave conditions corresponding 
to those of the Saco Bay wave-pattern photographs, refraction cases were computed using the 
three approximation techniques under investigation. For these cases, wave-ray and wave-crest 
patterns, as well as selected contours of equal wave height (constructed using a technique 
described in the appendix by James F. Kibler) are presented and compared. Selection of the 
appropriate techniques for bottom-topography approximation in other refraction • studies is 
discussed. 


SYMBOLS 

^ vector of constraints for constrained bicubic polynomial interpolation technique 

cjj coefficient of xV in approximating polynomial 

c vector of coefficients 

d depth, meters 

F matrix of independent variables evaluated at nodes of 16-node subset in 

constrained bicubic interpolation technique 

H wave height, meters 

h spacing between nodes in square grid 

P(x,y) bicubic interpolating polynomial 
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X matrix of independent variables evaluated at nodes of 12-node subset in 

least-squares techniques 

X spatial coordinate 

y spatial coordinate orthogonal to x 

Z general function of x and y 

e least-squares approximation error, meters 

X wavelength, meters 

Subscripts: 

k node index in x-direction 

C node index in y-direction 

m local node index in x-direction in 16-node subset 

n local node index in y-direction in 16-node subset 

o initial 

p pth row of X matrix 

A caret (^) denotes an estimate; a prime (') denotes a transpose of a matrix; and 
an underlined symbol denotes a vector. 

APPROXIMATION OF BOTTOM TOPOGRAPHY 

In general, computerized wave-refraction models employ bathymetry data which are 
supplied in a square-grid format with a fixed individual cell size. Used in conjunction with 
such a grid system are various techniques for approximating the continuous bottom- 
topography field to which the discrete bathymetry data apply. 
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Quadratic Least-Squares Approximation 

The model of reference 1 adopted a quadratic least-squares approach for topography 
approximation within the central cell of a localized 12-node subset of the total bathymetry 
field. This 12-node subset (see fig. 1) has also been shown to be the near-optimal con- 
figuration for use with a two-dimensional quadratic least-squares approximation (ref. 5). The 
approximating polynomial for the node in the subset located, for example, at coordinates 
(xjj,yg) can be written as 

=y^ x; + ^k,£ 

j=0 

where dj^ g is the actual depth at coordinates (xjj^,yg), cy are the six coefficients to be 
determined, and ej^ g is the approximation error at the node. For the entire subset, there 
are 12 linear equations of the form of equation (1) which can be written in matrix form as 

d = Xc + e (2) 

where d, c, and e_ are appropriate vectors and X is a matrix whose pth row contains 
values of the independent variables (xV^) evaluated at the pth node (p = 1,2, . . ., 12). 

The estimates of the coefficients resulting from application of least-squares principles 
can be written as 


c = (X'X)-^X'd (3) 

where c is the vector of estimates. Then the depth at any point (x,y) within the central 
cell of the 12-node subset can be approximated by 


2 2-i 


^(x,y) = X 2^ 
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CyxV 


i=0 j=0 


(4) 


Partial derivatives of the depth at point (x,y) can be approximated by evaluating analytical 
partial derivatives of equation (4) at that point. 

The matrix X (also X' and (X'X)*^) in this application contains terms which are 
functions only of the coordinates of the nodes of the subset of interest. Thus, by using 
local coordinates (that is, by measuring x and y with respect to some local origin 


5 


within the subset) the X matrix (and X', (X'X)‘^) can be made invariant over the 

entire region. With reference to equation (3), then, the only new information required for 
estimating the coefficients for any newly chosen 12-node subset is the vector of actual 
depths at those nodes. This mosaic least-squares approach can be applied very easily to a 
large geographical region. This approach also affords some smoothing of the input 
bathymetry data (which are prone to measurement error and, in reality, are subject to random 
geological perturbations), but offers no guarantee of continuity in the computed depths or 
derivatives between cells. 


Cubic Least-Squares Approximation 

Another technique which can be used for bottom-topography approximation is the 
cubic least-squares technique such as that employed in the model of reference 2. The 
approximating polynomial for the node in the bathymetry field located at coordinates 
(x]j,yg) can be written in this case as 

3 3-i 
i=0 j=0 

The model of reference 2 applied such a polynomial to the central cell of the general 
12-node subset shown in figure 1. For this subset a set of 12 linear equations of the form 
of equation (5) can be written similar to equation (2) (with appropriately changed 
dimensions). Least-squares estimates of the coefficients can be determined by an equation 
analogous to equation (3) (with 10 estimated coefficients rather than 6, as was the case with 
the quadratic approach). 

The depth at any point (x,y) within the central cell can then be approximated by 
3 3-i 

d(x,y) = CyxV-* (6) 

i=0 j=0 

Again, partial derivatives of the depth at the point (x,y) can be approximated by evaluating 
analytical partial derivatives of equation (6) at that point. 

Again the X, X', and (X'X)'^ matrices can be made invariant over the region 
through the use of local coordinates. Some smoothing of the input bathymetry data results, 
but to a lesser degree than the smoothing afforded by use of the quadratic least-squares 
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technique. In using the cubic least-squares method in a mosaic fashion, there again is no 
guarantee of continuity in computed depths or derivatives between cells. 

Constrained Bicubic Polynomial Interpolation 

The constrained bicubic polynomial interpolation technique has been used by the 
Coastal Engineering Laboratory of the U.S. Army Corps of Engineers for bottom-topography 
approximation. This procedure has the advantage of ensuring continuity in the depth and in 
certain derivatives of the depth over the entire region under study; its primary disadvantage 
lies in the fact that it does not provide any smoothing of the input data. However, a 
separate smoothing algorithm could be used to smooth the raw data prior to the use of the 
interpolation algorithm. 

The general form of the equation representing the approximating surface can be 
written as 


3 3 

P(x,y) = ^ 

i=0 j=0 

This polynomial is applied to the central cell of the 16-node subset shown in figure 2. The 
16 coefficients of the polynomial are determined by applying constraints on the polynomial 
and on certain' of its derivatives at the comer nodes of the central cell. The constraints 
also insure continuity in the polynomial and constrained derivatives between adjacent cells 
of interest. 

Four constraining conditions are that the value of the polynomial must equal the 
known depth values at the corner nodes of the central cell, or that 

^(^m’^n) “ *^m,n (For all combinations of m = 2,3; (8) 

n = 2,3) 


where m,n are local indices as shown in figure 2. 

Eight additional conditions can be stated by constraining the analytical first partial 
derivatives of equation (7) 9P/9x and 9P/9y at the corner nodes of the central cell to 

be equal to the central finite-difference approximations of these derivatives. These conditions 
can be written as 


^1 

x=Xjn,y=y„ 


^ %+l,n ~ ^m-l,n 


2h 


(For all combinations of m = 2,3; (9) 

n = 2,3) 
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and 


9y 




x=Xm,y=yn 


2h 


(For all combinations of m = 2,3; (10) 

n = 2,3) 


where h is the distance between adjacent nodal points. 

Four additional constraining conditions must be chosen in order to obtain a solution 
for the 16 polynomial coefficients. One could select from constraints involving any of 
the three second partial derivatives of the bicubic polynomial. However, constraining the 
mixed second partial derivative 0^p/dx 3y apparently affords some degree of smoothness 
in both coordinate directions along the approximating surface. Constraining either of the 
pure second partial derivatives 9^p/3x^ or 3^p/3y^ would guarantee smoothness in only 
one direction along the surface. Thus, the four remaining conditions are chosen so as to 
constrain the analytical mixed second partial derivative of equation (7), that is, 

3^p/3x 3y, to be equal to its finite-difference approximation at the comer nodes of the 
central cell. These conditions can be written as 


3^P 
3x 3y 


x=Xj„,y=yn 


‘^m+l,n+l ‘ ‘^m+l,n-l ' ^m-l,n+l ^m-l,n-l 

a 


( 11 ) 


4h^ 


(For all combinations of m = 2,3; 
n = 2,3) 


Conditions in equations (8) to (11), the left-hand sides of which involve forms of the basic 
polynomial (eq. (7)), thus constitute a set of 16 linear equations in the 16 unknowns 
cqO) cjq, . . ., C 23, C33. In matrix form, these can be written as 

Fc = a (12) 

where F is a 16 X 16 matrix with terms of the form xV (iJ = 0,1, 2,3), ^ is a 
16-element vector of the unknown coefficients cqq, Cjq, . . ., C23, C33, and ^ is a 
16-element vector the terms of which are the right-hand sides of conditions in equations (8) 
to (11). The solution to the matrix equation (12) is 

c= F'^a (13) 

The depth and its partial derivatives at any point (x,y) within the central cell can then be 
approximated by substituting these coefficients into equation (7) and into its appropriate 
analytical partial derivatives. 


8 



For a fixed grid system, the elements of the vector ^ are functions only of input 
depth values within the 16-node subset and the spacing between nodes. As in the least- 
squares solutions, if local coordinates are used, the matrix F is invariant as the bicubic 
interpolating technique is used mosaically over the region of interest. Since the 
matrix F is invariant over the region, F"^ is also invariant and must be computed only 
once for application to the region under study. Thus, solution for the coefficient vector 
for any 16-node grouping requires the supply of only the a vector, the terms of which 
again are functions only of input depth values at nodal points of the subset and the spacing 
between nodes. 


COMPARISON OF APPROXIMATION TECHNIQUES 

Since the three techniques for bottom-topography approximation discussed in the 
preceding section differ with regard to input data smoothing and to continuity of the 
approximating polynomial, quite different results might be expected to arise from the use 
of the different techniques in similar refraction calculations. In order to assess the magni- 
tude of differences in results, wave patterns and parameters, which are computed by using 
the different approximation techniques in the same refraction program with a fixed 
bathymetry data field and a fixed set of initial wave conditions, can be compared. 

For the present study, the basic computer routines of the mid-Atlantic continental- 
shelf wave-refraction program of reference 3 were adopted. As this basic program used 
only a quadratic least-squares topography approximation technique identical to that of 
reference 1, options were inserted to allow for the use of the cubic least-squares and 
bicubic polynomial interpolation techniques. In addition, the bathymetry data field for 
Saco Bay, Maine, was inserted in a 0.125-nautical-mile square-grid format in the program 
in place of the basic mid-Atlantic continental- shelf data. The Saco Bay region was 
selected for the present study for several reasons. The area is sufficiently small in size 
(9 n. mi. by 7.5 n. mi.) to allow for much easier comparison of overall wave patterns and 
distributions of parameters such as wavelength and wave height than would be possible for 
larger regions. At the same time, the complex nature of the Saco Bay bottom topography 
and shoreline leads to complex wave patterns which are representative of patterns experi- 
enced in larger near-coastal regions. A three-dimensional computer-drawn plot of the 
Saco Bay bathymetry data field (with a vertical exaggeration of approximately 75 to 1) is 
shown in figure 3. In addition, computed wave patterns can be referenced to an actual 
wave pattern which was observed in vertical aerial photographs of Saco Bay (fig. 4) and 
which was used to provide initial conditions and a qualitative verification of the refraction 
calculations. 
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RESULTS AND DISCUSSION 


In order to compare the topography approximation techniques discussed previously, 
three cases were computed using the modified mid-Atlantic continental -shelf refraction 
program with the Saco Bay bathymetry field. From a photogrammetric analysis (ref. 6) of 
the wave pattern seen in figure 4, the initial wavelength Xq to be used in each of the 
cases was determined to be approximately 100 meters; the initial direction of wave propa- 
gation was fixed at -135° (with respect to the positive X-axis). By using linear wave theory 
with an assumption of deep water, the initial wave period was determined from the 
100-meter initial wavelength to be 8 seconds. Since the initial wave height could not be 
determined from the photograph and since it is not a critical parameter with respect to 
computation of the wave pattern itself, a reference initial height of 1.22 meters was 

arbitrarily selected. 


Wave Patterns 

Computed wave-ray (orthogonal to wave crest) patterns for the three cases are 
presented in figure 5. In a comparison of the three computed patterns, obvious differences 
in the paths of travel of individual rays can be seen. Of particular interest is the path of 
ray number 36 in the pattern computed using bicubic interpolation. This ray, which is 
highlighted in figure 5(c), exhibits an abrupt change in direction at coordinates x = 15 grid 
units and y = 14 grid units, near which point (see fig. 3) is located an isolated extremum 
in the bathymetry data. This abrupt change in direction illustrates the sensitivity of the 
refraction calculations to isolated extrema (whether real or erroneous) in the bathymetry 
data. In turn, such extrema in the bathymetry field are emphasized with the use of an 
interpolative technique such as the present constrained bicubic polynomial technique. This 
sensitivity suggests that careful verification of real extrema and the removal (such as by 
presmoothing the data) of erroneous extrema from the data field must be done in order to 
avoid possibly misleading computations with the use of the constrained bicubic interpolation 
technique. Neither of the patterns computed using the least-squares theory (figs. 5(a) and 
5(b)) exhibits such abrupt changes in ray directions, since the least-squares approach tends 
automatically to smooth isolated disturbances in the data field. 

However, from an overall viewpoint the three patterns show a great deal of similarity. 
On aU three patterns, rays in the number range from 15 to 30 tended to converge and to 
stop near the island located near x = 26, y = 36. This tendency is indicative of waves 
breaking as they approach the island. (See fig. 4.) Similarities also can be noted among 
the three patterns in the groups of rays which refract around the island, although the group 
of rays numbered 30 to 40 is less diffuse in the quadratic least-squares pattern (fig. 5(a)) 
than in either of the other two patterns. 
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Wave-crest patterns can be derived from computed ray patterns by connecting points 
which occur at the same time on adjacent rays. The points must be connected in such a 
manner that the connecting line (crest) is orthogonal to each ray. Accordingly, ray patterns 
identical to those of figure 5, with the exception of tic marks inserted to denote equal 
time, were computed. Wave-crest patterns derived from these ray patterns are shown in 
figure 6. A great deal of similarity can be seen among the three patterns, as could be 
expected considering the similarities in the ray patterns of figure 5. The most prominent 
features in the crest patterns are the refraction of crests around the island located 
near x = 26, y = 36 and the interaction of two distinct trains of refracted crests in the 
region near x = 15, y = 25 to 30. A qualitative examination of the Saco Bay vertical 
aerial wave photographs (fig. 4) shows similar features in the actual wave pattern. 

Contours of Wave Parameters 

A more quantitative assessment of the effects of using different bottom-topography 
approximation techniques in refraction calculations can be made by comparing spatial 
distributions of computed wave parameters. Two of the more fundamental parameters 
which can be compared in this manner are the wavelength X and the wave height H. 

An effective means for presenting spatial distributions of such parameters is the construction 
of contours of equal values for the parameter of interest. 

Contours of selected values of wavelength and wave height were constructed for the 
three cases under study by using a technique described in the appendix. The arrays of 
data points supplied to the contouring program included the wavelength and wave height 
(and the spatial coordinates) at every odd-numbered data point (in time) along each of the 
40 computed rays. Contour curves were then constructed using the cubic spline curve- 
fitting technique with a tension factor of -3. (See appendix.) 

Contours of wavelength for levels of 25, 50, and 75 meters are presented for the 
three topography approximation techniques in figure 7. Similarities between distinct groups 
of contours can be noted among the three regional plots. One such group comprises the 
closed contours for wavelength levels of 75 and 50 meters located near x = 25, y = 36. 
Smaller closed contours at the 25-meter level lie interior to these 50- and 75-meter 
contours on the plots corresponding to the cubic least-squares and constrained bicubic 
interpolation techniques (figs. 7(b) and 7(c)). These grouped contours are indicative of the 

reduction in the length of the wave as it approaches and refracts around the island located 

near these same coordinates. The other distinct group comprises the curved open contours 
at all three wavelength levels which lie nearer the shore and points out the gradual 
reduction in the length of the wave as it travels over increasingly shallower water in 

approaching the shore. As a result of the more constricted local ray pattern (fig. 5(a)) 

computed using the quadratic least-squares method, the open contours corresponding to the 
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quadratic least-squares approximation technique (fig. 7(a)) do not approach the X-axis as 
closely as do those contours corresponding to the other techniques. 

The computed wave height is most strongly affected by the relative separation of 
individual wave rays passing through the local area of interest. (See ref. 7 for details of 
computation.) For example, as adjacent rays converge, the wave height at points between 
the two rays tends to increase, and vice versa. Since differences were seen in the paths of 
corresponding rays in the patterns computed using the different topography approximation 
techniques, local differences would be expected in comparative wave-height contours. For 
example, in comparing height contours at the 2.0-meter level for the different techniques 
(fig. 8), local differences can be noted. However, the contours in the quadratic and cubic 
least-squares diagrams tend to lie in a pattern oriented in the initial direction of motion of 
the rays and located on either side of the islands. No such definitive pattern is seen in 
the bicubic interpolation diagram, although the contours are located seaward and to either 
side of the island. Similar trends can be noted in the contour diagrams at the 1.5- and 
1.0-meter levels which are shown in figures 9 and 10, respectively. However, the pattern 
in the cubic least-squares diagrams (figs. 9(b) and 10(b)) tends to deviate somewhat from 
the nominal pattern with the appearance of additional contours shoreward of the island 
near the center of the general pattern. Comparative contours at the 0. 5-meter wave-height 
level are presented for the three approximation techniques in figure 11. Although great 
differences can be seen among the three diagrams, the contours are grouped shoreward of 
the island and, in several areas, are oriented in the same general direction. 

Overall Comparison 

In an overall manner, the computed wave-ray and wave-crest patterns are not strongly 
affected by the choice of different topography approximation techniques. Comparative 
patterns agree equally well in a qualitative manner with vertical aerial photographs of the 
reference Saco Bay wave pattern from which initial conditions for the comparative cases 
were drawn. Comparative contours of wavelength and wave height were constructed using 
computed data along individual wave rays and these contours exhibit similar characteristics 
in overall contour orientation. However, comparative height contours, in particular, show 
considerable local differences. In the absence of more quantitative data on the spatial 
distribution of wavelength and wave height, the best approximation technique cannot be 
defined; the choice of any of the three techniques could be made with a similar degree of 
confidence in the resultant overall computed patterns. For those applications in which local 
accuracy in computed results is critical, more extensive quantitative comparison of computed 
results with wave-pattern photographs and wave-height data (as might be obtained from use 
of radar or a laser profilometer) would be required in order to define the best approxima- 
tion technique. 
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CONCLUDING REMARKS 


A study of the effects of using different techniques for approximating bottom 
topography in a wave-refraction computer model has been conducted. The basic Lan^ey 
Research Center — Virginia Institute of Marine Science wave-refraction model was modified 
to study the region of Saco Bay, Maine, and to include optional topography approximation 
techniques involving a quadratic least-squares approach, a cubic least-squares approach, or a 
constrained bicubic polynomial interpolation approach. Comparative refraction cases were 
run using the three approximation techniques with ah assumed initial wave height and fixed 
initial conditions for wave period and direction derived from vertical aerial photographs of a 
Saco Bay wave pattern. 

Wave-ray and wave-crest patterns computed using the three techniques were quite 
similar from an overall standpoint, and the crest patterns compared well, in a qualitative 
sense, with the reference photographs. However, some differences were noted in the paths 
of corresponding rays in the comparative patterns. The path of one particular ray, which 
was computed using the constrained bicubic interpolation technique, exhibited an abrupt 
change in direction in the vicinity of an isolated extremum in the bathymetry data field. 
This change in direction illustrates the sensitivity of refraction calculations to such an 
extremum (whether real or erroneous) the importance of which, in turn, is emphasized with 
the use of an interpolation technique such as the present constrained bicubic polynomial. 
This sensitivity implies that verification of real extrema and removal of erroneous extrema 
in the data field are required in order to avoid possibly misleading computations with the 
use of the bicubic interpolation technique. Contours were constructed at selected levels of 
wavelength and wave height for the three comparative cases. Comparative wavelength 
contours were quite simUar, and even though considerable differences could be observed in 
local areas in the comparative wave-height contour diagrams, a similar pattern of orientation 
on a gross scale could be noted. 

In light of the absence of actual wave-height data for comparison, the best approxima- 
tion technique could not be defined, and any of the three techniques could be expected to 
yield comparable overall results. For applications in which local accuracy in computed 
results is critical, more extensive quantitative comparison of computed results with actual 
wave data would be necessary in order to determine the proper topography approximation 
technique. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, Va. 23665 
August 14, 1975 
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APPENDIX 


DESCRIPTION OF CONTOURING TECHNIQUE 

James F. Kibler 
Langley Research Center 


A contour plot is a convenient way to display a function of two independent 
variables. The contour is simply a line connecting constant values of the function plotted 
on a grid composed of the corresponding independent variables. In this way, a three- 
dimensional surface can be easily visualized in two dimensions. The purpose of this 
appendix is to describe briefly the contouring technique used to generate the contour plots 
of wavelength and wave height presented in this paper. The technique was developed by 
C. L. Lawson of the Jet Propulsion Laboratory and modified by A. K. Cline of the 
Institute for Computer Applications in Science and Engineering at the Langley Research 
Center. 

Assume that the function to be contoured is available as an array of randomly 
spaced data points. Each data point consists of an x, y, and Z value. The x and 
y values correspond to the independent coordinates to be plotted and Z = Z(x,y) is 
the function value to be contoured. The contouring process begins by forming a convex 
region composed of triangles from the set of x-y data points (fig. 12). Figure 12(b) 
illustrates an arrangement of triangles which form a convex region from the data points in 
figure 12(a). Obviously, this arrangement is not unique. A somewhat arbitrary criterion 
is established to maximize the smallest interior angle between possible pairs of triangles. 

For example, with a set of four data points there are two choices for arranging the two 
triangles (fig. 13). Under the above criterion, figure 13(b) is chosen as the preferred 
arrangement. That is, the smallest interior angle in figure 13(b) is larger than the smallest 
interior angle in figure 13(a). 

After the triangles are arranged, each Z value is assigned to its corresponding x-y 
vertex- in each triangle. An interpolating plane which fits the vertices exactly is then 
calculated for each triangle. If a chosen contour level Ues between the Z values for 
any triangle, the two contour points which lie on the sides of the triangle are computed 
from the interpolating plane (fig. 14). These contour points are calculated for each 
triangle in the convex region and for each specified contour level. The contour points for 
any given contour level may then be connected by straight lines to yield a polygonal 
approximation to the required contour. Alternatively, a method involving a cubic spline 
function under tension can be used to fit a smooth curve to each set of contour points. 
(For example, see ref. 8.) This curve fits the interpolated contour points exactly and 
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varies smoothly between them. The tension factor may be adjusted to yield suitable 
curves. (A large factor corresponds to straight lines between contour points.) 

The contours which result from this technique are only approximate. Their accuracy 
is limited by the linear interpolating procedure and is dependent upon the density of the 
available data. 
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Figure 1.- Twelve-node subset of bathymetry data field to which quadratic and cubic 

least-squares approximations are applied. 
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Figure 2.- Sixteen-node subset of bathymetry data field to which constrained bicubic 

polynomial interpolation is applied. 
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Figure 4.- Mosaic of vertical aerial 
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(a) Quadratic least squares. 

Figure 5.- Saco Bay wave-ray patterns computed by using different topography 

approximation techniques. 
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(b) Cubic least squares. 
Figure 7.- Continued. 
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(c) Constrained bicubic interpolation. 
Figure 8.- Concluded. 
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(b) Cubic least squares. 
Figure 11.- Continued. 
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(a) A set of input data points. 
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(b) A convex region composed of triangles. 
Figure 12.- Contouring process. 
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Figure 13.- Two choices for arranging triangles for four data points. 


43 




Figure 14.- A triangle with linearly interpolated contour points. 
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